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PREFACE 


Since the problem of two bodies is the only problem in astrodynam- 
ics with a known solution for aribtrary initial conditions, it has been used 
in an approximate solution to the restricted problem of three bodies in the 
form of patched conic orbits. Since the development of the patched conic 
technique, several methods of approximating the solution to the restricted 
problem of three bodies have been presented, but none of them utilize full 
knowledge of the known integrals for the exact motion. It is believed that 
a method that uses knowledge of the known functions of the motion and is 
conceptually simple would be quite useful for studies of future space missions 

This study presents a method of calculating trajectories for the 
restricted problem of three bodies using conic motion that is frequently 
corrected in position and velocity. The correction in position and velocity 
is calculated using knowledge of the existing integrals or slowly-varying 
functions of the motion. This method is easily described. Assume that the 
trajectory has just been corrected. The motion to the next correction point 
and the correction there will be described. The independent variable is the 
magnitude of the radius vector. A change in the independent variable Ar 
is chosen and the trajectory is conically advanced through the interval Ar . 
Since the value of the function of the motion evaluated on the conic trajec- 
tory is not the same as the value predicted for the exact motion, position 
and velocity corrections are applied to the conic trajectory so that the 
value of the function will be the same as the predicted value. The process 
is repeated until the terminal conditions are reached. 

The results of this method are compared with numerically integrated 
trajectories. This method is qualitatively compared with other methods of 
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solution for the restricted problem of three bodies. 
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ABSTRACT 


This report presents a method of calculating trajectories for the 
restricted problem of three bodies which utilizes conic propagation of the 
state vector with frequent correction of position and velocity by means of 
a constant or slowly-varying function. This fast and accurate method of cal- 
culating trajectories has been applied to the planar circular restricted 
problem of three bodies , the planar elliptic restricted problem of three 
bodies 5 and the ephemeral restricted problem of three bodies. Two methods 
(the "refined” method and the "straight-forward" method) of determining the 
direction of the position correction (n^) are presented for the circular 
restricted problem and the elliptic restricted problem of three bodies. Only 
the "straight-forward" method is used with the ephemeral restricted problem 
of three bodies. The Earth, the Moon and a space vehicle comprise the res- 
tricted three body model that is used. Earth-to-Moon trajectories with per- 
ilune altitudes varying from 59 to 4551 nautical miles are calculated and 
compared at perilune with numerically integrated and patched conic trajec- 
tories. The results, as compared to the numerically integrated trajectories, 
are within 0 . 2 % in position and velocity vector magnitude (relative to the 
Moon) for the "straight-forward" and the "refined" choices of the position 

correction direction (n ) . 

c 

A detailed discussion of the two methods of choosing n c is pre- 
sented. A qualitative comparison between this method and other methods of 
calculating trajectories for the restricted problem of three bodies is also 
presented. 


IV 
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NOMENCLATURE 


The following list presents all significant symbols and abbrevia 
tions used in the main body of the text. Each symbol is accompanied by a 
brief description and the number of the equation where the symbol is intro 
duced. 


Vectors : 
A 



n 

c 


n 

P 


R, 





r 


1 


_ o 

vector (A = r) used in Equation (2.36) 
average perturbing acceleration (2,33a) 
vector used in Equation (2.37) 

unit vectors relative to the space vehicle (3.7) 


unit vectors describing primary two relative to primary one 
(2.13) 

angular momentum of massless particle (2.18) 

unit vector in position correction direction (2.29) 

unit vector in the velocity correction direction (2.23) 

position vector of primary one relative to CM (2.3) 

position vector of primary two relative to CM (2.3) 

position vector of primary two relative to primary one 

position vector of massless particle relative to CM (2.1) 
position vector of massless particle relative to primary one 
( 2 . 1 ) 


x 



r 2 position vector of massless particle relative to primary two 

( 2 . 1 ) 

AR change of position vector (7.2) 

AV change of velocity vector (7.1) 

6r position vector correction (2.28) 

<5r velocity vector correction (2.28) 

w angular velocity of primaries about CM (2.10) 


Scalars : 


E 


2/1 


e 


2/1 


2/1 


G 


J 

J 



m 


2 


R 


1 


R 


2 


r 


r 


1 



eccentric anomaly of primary two relative to primary one (4.2) 

eccentricity of primary two relative to primary one (4.2) 

angular velocity of primary two relative to primary one (4.3) 

universal gravitational constant in Section (3.4) 

Jacobi function (2.16) 

rate of change of Jacobi function (2.17) 

J evaluated on the conic trajectory (2.27) 

dimensional mass of primary one (2.2) 

dimensional mass of primary two (2.2) 

projection of r^ onto the Earth-Moon plane for the ephemeral 
restricted problem of three bodies 

projection of r^ onto the Earth-Moon plane (5.7) 

magnitude |r| (2.1) 
magnitude | r | (2.1) 

magnitude |r 2 | (2.1) 


xi 



v initial value of r (2.25) 

r~ final value of r (2.25) 

^ i 

■i 

Ar increment in r (2.25) 

Ar initial increment in Ar (2.25) 

o 

Ar_^ final increment in Ar (2.25) 

At increment in time (t^ - t^) (2.26a) 

6t time correction (2.30) 

<$r magnitude |6r| (2.29) 

6r magnitude, 1 6 r | (2.29) 

a angle from e to e [Figure (3.2)] 

X IT* 

* mrr 

0^ (i = 1,2) angular velocity of r^ 

y mass ratio-parameter (2.2) 

a variable = ±1 (3.11) 

d> angle between e and n (3.7) 

■ , r c 

to magnitude | to] 


Subscripts : 

f indicates final value 

i indicates 1 or 2 

o indicates initial value 

p indicates perturbing acceleration direction 

x referenced to x direction 

y referenced to y direction 

z referenced to z direction 
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Abbreviations : Numbers to the side represent Sections where they first appear. 

CM center of mass (2.2) 

er Earth radii (3.4) 

fps feet/second (3.4) 

hr hour (3.4) 

min minute 

n. mi nautical mile (3.4) 

rad radians (Table 1) 

sec second (3.4) 
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CHAPTER 1 


INTRODUCTION 

1.1 General Background 

Since the problem of two bodies is the only problem in astrodynamics 
with a known solution for arbitrary initial conditions , it has been used ex- 
tensively as a model of the problem of planet -orbiting satellites. The solu- 
tion to the problem of two bodies is also a good approximation to the motion 
of the planets relative to the sun. However, it is not a good approximation 
for Earth-Moon or interplanetary trajectories, because it cannot include mul- 
tiple force centers . This led to the use of the restricted* problem of three 
bodies as a mathematical model for Earth-Moon trajectories and successive 
portions of interplanetary trajectories. 


1.2 Restricted Problem of Three Bodies 

The restricted problem of three bodies in this study is the motion 
of a massless particle (space vehicle) in the vicinity of two massive primar- 
ies (see Figure 1.1). The’ unit base vectors e^ and e are in the plane 



Figure 1.1 Diagram for the Restricted Problem of Three Bodies 


^Restricted in the sense that the mass of the third body is small enough that 
it does not affect the motion of the two primaries. 


1 



2 


containing the motion of the two primaries relative to the center of mass 

(CM) and e is the unit vector perpendicular to that plane, 
z 

1.2.1 Planar Circular Restricted Problem of Three Bodies 

The circular restricted problem of three bodies is the configuration 
where the two primaries are in circular orbits about the CM and the motion of 
the particle is in the plane of motion of the two primaries. 

1.2.2 Planar Elliptic Restricted Problem of Three Bodies 

The elliptic restricted problem of three bodies is the system where 
the two primaries are in elliptic orbits relative to the CM. The motion of 
the particle is again in the plane of the motion of the two primaries. 

1.2.3 Ephemeral Restricted Problem of Three Bodies 

The ephemeral restricted problem of three bodies is the three-dimen- 
sional motion of the massless particle, where the position and velocities of 
the primaries (the Earth and the Moon) are obtained from available ephemeris 

information. Such ephemeris information is available in readily accessible 

TnT , . [6] * 

form for computer use on the JPL ephemeris tape 

For arbitrary initial conditions there are no known analytic solu- 
tions to any of the above mentioned problems. Since the restricted problem 
of three bodies is a representative mathematical model of the Earth-Moon space 
vehicle system and of successive parts of interplanetary trajectories, it is 
desirable to have a fast, accurate solution from the standpoint of guidance 
and trajectory analyses. This solution can be used to determine parameter 
sensitivity and guidance sensitivity for several trajectories with little 

^Numbers appearing in the text as superscripts indicate references listed in 
the Bibliography. 
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computer time expense. To perform a similar analysis using a numerical inte- 
gration technique would be very expensive in terms of computer time. 

1.3 Approximate Solutions 

[7] 

The patched conic, introduced by Egorov in 1958, was one of the 
first approximate solutions to the restricted problem of three bodies. The 
patched conic for the restricted problem of three bodies consists of two 
conic segments, the conic of a particle about primary one without the pertur- 
bations of primary two and the conic of the same particle about primary two 
without primary one perturbations, which are joined at a point in space to 
produce the composite trajectory. The joining point in space is taken to lie 
on the surface of a nearly spherical surface , centered at primary two , which 
is called the Mean Surface of Influence and is discussed in Ref. [8] (see 
Figure 1.2). 


Mean Surface of Influence 



Hyperbolic 
Segment 
relative 
to Primary 


£ 


Figure 1.2 Patched Conic Geometry for Primary 1 to Primary 2 Trajectories 
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Some of the disadvantages of the patched-conic are: 

1. it produces large errors for trajectories that have long transit 
times 3 

2. it is sensitive to the choice of the magnitude of the mean sur- 
face of influence , and 

3. it has no means of including the effect of the perturbing body. 
Its advantages are : 

1. it is very fast computationally 

2. it provides reasonable velocity requirements if the initial and 
final position vectors are given. 

[11 12 13] 

The method of matched asymptotic expansions 5 5 is another 

method of approximating the trajectory of a particle in the presence of two 
primaries. The initial difficulty of the matched asymptotic expansion is 
the algebra and computer program check out required to obtain suitable re- 

[14] 

suits. Until later refinements were applied to the method , the computer 

time required to obtain a solution is almost as large as the time required to 

obtain a numerically integrated trajectory. 

During the period of time from 1963 to 1966 much work was done at 

NASA and TRW Systems to improve the patched-conic by applying a velocity cor- 

[27 28 ] 

rection at the patch point ’ . The velocity correction being calculated 

from knowledge of the "Jacobian Function" for the restricted problem of three 
bodies. This method was an improvement to the patched conic, but was not suf- 
ficient for all cases. 

During the period of time from 1967 to 1969 the Hybrid Patched-Conic 

fg 1 

Technique was developed by Escobal, et_ al_ . At first appearance it seemed 

that the Hybrid Patched Conic Technique was accurate and fast enough to meet 
the needs of NASA and industry at the time. The disadvantages to the Hybrid 



5 


Patched Conic Technique are that it requires a patched conic solution for a 
reference trajectory, and its accuracy is limited if the perilune altitude 
is large (e.g., greater than 3000 n.mi.). 

In 1969 this investigation was initiated using knowledge of the 
"Jacobian Condition” and the "angular Momentum condition" to make corrections 
at several points along the trajectory to see if this would not produce a 
quickly-calculated trajectory that was sufficiently accurate, as compared to 
a numerically integrated solution. . It was later determined that the angular 
momentum correction was not accurate enough to help improve the accuracy. 

After the present investigation was initiated, it was learned that 
several individuals at TRW Systems at Houston were working on the same type 
of problems but with quite different approaches. They developed the multi- 
conic method^ 4 ^ and the pseudo-conic method^ 

The multi-conic method uses two-body motion as the basic propaga- 
tion technique. Gravitational effects are accounted for by assuming that 
each perturbing body causes' independent two-body motion. The effects are 
then summed along the trajectory. The procedure does involve a retracing 
step and using a zero gravity step. Thus, the method is more, complicated 
than the "Jacobian" correction method presented here. 

The pseudo-conic method also uses conic motion as method of propa- 
gation, but it continues on past the mean surface of influence along a tra- 
jectory that is regarded as a pseudostate . Then it propagates from the mean 
surface of Influence to the desired final time. The pseudo-conic does reduce 
the patched conic error considerably, but it does not seem to be as accurate 
as the multi-conic method. 

An "integral hypersurface" technique^ 1 ^as ^ een use< ^ by 

Tl7l 

Nacozy to constrain the numerically integrated solution to remain on the 
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integral surfaces. A similar technique has been used iteratively by Miller^ 15 "^ 

in a gravitational n-body integration to control the usual ten first integrals 

rigi 

of motion. Miller also used the first ten integrals of the equations of 
motion as controls for nobody integration. In a comparison of a corrected 
solution of the system with a similar , uncorrected solution, he finds that 
the two solutions diverge from each other - indicating the instability of the 
gravitational system. Aarseth^^ used a similar integral surfaces technique 
to correct the integrals, the positions, and the velocities of the computed 
solution to account for the removal of escaping bodies from the system. 

1.4 Motivation 

The motivation for the approach taken here is that proper use of 
the knowledge obtained from the "Jacobian" Function could produce results that 
are a significant improvement over the results obtained from a patched conic 

trajectory with much less computer time than is required for numerical inte- 

. . [24] 

gr at ion 

The first step was to apply the theory to the planar circular res- 

[23] 

tricted problem of three bodies. Since the Jacobian Function is a constant 

for the circular restricted problem of three bodies , it was felt that it would 
be best to apply the theory to the circular restricted problem of three bodies 
before proceeding to the elliptic and ephemeral restricted problem of three 
bodies. The description of the method and the necessary equations are de- 
rived in Chapter 2. The application to and the results obtained from the cir- 
cular restricted problem of three bodies are presented in Chapter 3. Next, 
the elliptic restricted problem of three bodies is treated and the results 
presented in Chapter 4. Last, the ephemeral restricted problem of three bodies 
application and results are presented in Chapter 5. A detailed discussion of 
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the choice of the position vector correction direction (n ) and the velocity 
vector correction direction (n^) P resen "^ e< ^ i n Chapter 6. The summary 

and conclusions are presented in Chapter 7. 

A qualitative comparison of the patched conic, the hybrid patched 

conic technique, the matched asymptotic expansion technique, the multi-conic,. 

i 

and the pseudo-conic is also presented in Chapter 7 . 



CHAPTER 2 


DESCRIPTION OF METHOD 


2.1 Assumptions 


If the effects of the gravitational fields of the sun and other 
planets are neglected, the system containing the Earth, the Moon, and a space 


vehicle can be approximated as a three-body system (see Figure 2.1), 




MOON 


Figure 2.1 Configuration of the Earth, the Moon, and the Space Vehicle. 

This can be modeled, as in Figure 2.2, as the restricted problem of 
three bodies. The restriction is that the space vehicle (or the massless 
particle) does not affect the motion of the two primaries. The masses of the 
two primaries are assumed to be spherically symmetric and homogeneous in 
concentric layers. 


8 
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e y 



Figure 2-2 Diagram for the Restricted Three-Body Problem (R3BP). 


The problem will be formulated in a vectoral notation that can be 
used for both two- and three-dimensional problems. For the ephemeral res- 
tricted problem of three bodies the position and velocities of the primaries 

r c "I 

(the Earth and the Moon) are obtained from the JPL ephemeris tape L 


2.2 Development of the General Equations of Motion 

The non-dimensional equation of motion for the massless particle is 


.. r r 

r+(l-M)— + y — | = 0 


( 2 . 1 ) 


(see Figure 2.2) where y , the mass ratio-parameter, is 


nu 


y = 


(m + m 2 ) 


( 2 . 2 ) 


and m^ and m^ are the dimensional masses of the two primaries. 

r = acceleration of the massless particle relative to the center 
of mass (CM) . 

r = position vector of the massless particle relative to the CM. 
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r^ = position vector of the massless particle relative to primary 


one . 


r^ = position vector of the massless particle relative to primary 
two. 


Therefore, 


= r - R, 


= r - R. 


, r 2 = r - r 2 


, r 2 = r - r 2 


( 2 . 3 ) 

( 2 . 4 ) 


R 1 " yR 2/l 


R 2 = (l - y)R 2/l 


R 2/l " R 2 " R 1 ®x R 2/l 


( 2 . 5 ) 

( 2 . 6 ) 
( 2 . 7 ) 


where 


R 


2/1 


R. 


= Position vector of Primary two relative to Primary one. 
= Position vector of Primary one relative to the CM. 


R^ = Position vector of Primary two relative to the CM. 

The equations of motion relative to primary one and primary two are, respec- 
tively, 


? 1 + (1 - y) + y i 

r i 


r- 

r_ 


R 


2/1 


3 + -3 
r 2 R 2/l 


= 0 


( 2 . 8 ) 


and 


r 2 + 4 — 
r 2 


+ (1 - li) 


1 _ 

3 


R 


2/1 


R 


2/1 


( 2 . 9 ) 


The velocity r and acceleration r can be written as 
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r = r + a) x r 


00 _ _ o •_ 

r + u) x (u) x r) + 2a)xr+o)xr 


(2.10) 

(2.11) 


where 


? - St® 


(2.12a) 


r =■ dr/dt with e , e , e fixed. 

x y z 


(2.12b) 


In the cartesian coordinate system indicated by the rotating, unit base vec- 
tors e , e , e , 
x y z 


r = xe + ye + ze 
x J ' y z 


0 

r 


xe + ye t ze 
x y z 


(2.13) 

(2.14) 


a) = angular velocity of primaries about the CM . 


= we 


(2,15) 


The Jacobian function for the restricted problem of three bodies and its 

, . .. [3,29] 

derivative are 

1 ° 2 1 _ __ 

J = — r • r - — (to x r) • (w x r) - (1 - u)/*^ - u/r 2 (2.16) 


and 


J = -co • h + y(l - y)R 


2/1 


(2.17) 


where 


= angular momentum of the massless 
particle relative to the CM . 

_ I _ 2 _ 9 - - - 

= rxr = rxr + tor - r (w • r) 


J> 


(2.18) 
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This assumes that no other forces are acting on the system and 

Ki + rr 1 = 0 


(2.19) 


R 


2/1 


which implies that 


w - u>e z ( R 2/l x R 2/l^ /R 2/l 


( 2 . 20 ) 


For the case of the circular restricted problem of three bodies 
it is clear that 


to = 0 and R, 


2/1 


0 


( 2 . 21 ) 


This then leads to the well known Jacobian integral for the circular res- 
tricted problem of three bodies (see pp. 16 of Ref. 29). That is, 
o o _ 

J = — r • r - — (to x r) • (to x r) - (1 -y)/r^ - y/r^ - Const . (2.22) 

Further details of each of these equations (2.1-2.22) will be dis- 
cussed as necessary in the remaining Chapters. 

It is desirable at this point to discuss the method of application. 


2.3 Method of Calculation of Trajectories 

Due to its computational simplicity, conic motion has been chosen 
to be the method of trajectory advancement. The force center is the primary 
on the same side of the surface of influence as the massless particle (the 
surface of influence is defined on p. 148 ff of Ref. 8) and the independent 
variable is the magnitude of the position vector from the force center. 

This choice of independent variable eliminates the need for iteration involv- 
ing Kepler's equation. 

After the trajectory has been conically advanced over the desired 
position vector magnitude interval, a correction to the position and velocity 


I 
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is calculated using the Jacobian function, which is constant or slowly-vary- 
ing for the exact motion. Slowly-varying functions [Equation (2,17)3 must be 
integrated over the propagation interval. 

In the application of this method to the restricted problem of three 
bodies one scalar (Jacobian) function is involved in the correction procedure. 
This function is used to correct one velocity vector component and one posi- 
tion vector component. The direction of the velocity component is the approx- 
imate direction of the time-averaged perturbing acceleration. The position 
component T s direction is different and is discussed in Chapter 6. 

The direction of the perturbing acceleration is obtained from Equa- 
tions (2.8) and (2.9). For motion relative to primary one the unit vector in 
the direction of the average perturbing acceleration 


n 

P 


- ( K 2 / 1 / R 2/1 




R 2/l /R 2/l “ V r ? 


(2.23) 


where | | indicates the absolute value. For motion relative to pri- 

mary two this is 


n 

P 



- R 2/l /R 2/l + V r !] 

R 2/l /R 2/l + f l /r l 


(2.24) 


Further details of the correction direction procedure will be given 
as necessary in the appropriate sections. 

Assume that the trajectory has just been corrected. The motion to 
the next correction point and the correction there will be described. 

The interval of propagation Ar is chosen to vary linearly with r 
and is calculated by means of the equation. 
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Ar = Ar + (At* - Ar )(r - r )/(r - r ) 

o f o o r o 


(2.25) 


Where r and Ar are initial values and rv and At* are final values, 
o o r r 

The state vector is conically propagated to a new state vector at r + Ar . 

Due to the choice of independent variable, no iterations are necessary and 

this is a straightforward procedure. 

At the new conically-advanced state, a "conic” Jacobian function 

and its derivative J and J are calculated. The approximate values of 

c c 

J and J at the newly advanced state are predicted using trapezoidal inte- 
gration. That is. 


= t 0.5(J^ + J^At 


(2.26) 


where J 2 is the predicted value of J at r + Ar , and are the 

values of J and J evaluated at the previous r , and 


At = t 2 - t x 


(2.26a) 


where t is the time associated with the trajectory at r and t g is the 

• • 

time associated with the trajectory at r + Ar . The derivative J 2 is J 
evaluated on the conic trajectory at t . Then J J 2 . Since the exact 


motion is not conic 


J i J 
c 


(2.27) 

At this point, correct both velocity and position vectors. 

r> -y r + Sr j r r + Sr (2.28) 


and 


Sr = n Sr ; Sr = n Sr 
p c 


(2.29) 


The unit vector n indicates the direction of the velocity correction while 
P 
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n indicates the position vector correction direction. These are discussed 
c 

in Chapter 6. 

Note that, since the time correction 

fit S 0 (2.30) 

the following equations are true: 

• • • • 

Sr = 5r^ = Sr 2 ; Sr = 6^ = Sr 2 (2.31) 

The equation for determining Sr and Sr is 

1 ° 0 £ £ i £ £ i _ 

— (r + 6r) ’ • (r + 6r) - — r ' r - - tw x (r f <$r)] 

* [5 x (r + Sr)] + ~ ( a) x r) • (u) x r) - (1 - y)/(r^ + (2.32) 

+ (1 - v)/r 1 - y/(r 2 + 6r 2 ) + y/r 2 = J “ J c 

Since two scalar quantities are being corrected by means of one scalar func- 
tion, a relationship between Sr and Sr is needed. If a is the average 

perturbing acceleration, 


6r 

= aAt 


(2.33a) 

Sr 

= i a(At) 2 = I ( Sr) (At ) 


(2.33b) 

where At is the time interval corresponding to Ar and is 
Kepler’s equation. Ignore the difference between n^ and 
6) and use 

evaluated from 
n^ (see Chapter 


Sr = j (Sr) (At ) 


(2.34) 

Equation (2.32) 

• 

is linearized with respect to Sr 

and Sr 

. The 

resulting expressions for 

• 

Sr and Sr are 
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with A 
B 


J - J 

r * c 

6r . = — I I J I I 

[(n • A) + |- (n • B ) ( At ) ] 

p 2 c 

- (J - J )At 

r 1 c 

6r = 2 I 1 )~Z I 

1 E(n • A) + i (n • B ) ( At) ] 

p 2 c 

o 

r 

_ 2 2--- -3-3 

{ (a> x r) - [a) - (a) * n )]r + (l - y)r /r + yr /r } 

C 11 z z 


( 2 .35a) 


(2.35b) 


(2.36) 

(2.37) 


After the corrections 6r and 6r are made to the trajectory , . the 
process is repeated until the desired stopping condition (perilune) is sat- 
isfied. 


The method is applied to the circular restricted problem of three 
bodies, the elliptic restricted problem of three bodies and the ephemeral res- 
tricted problem of three bodies in Chapters 3, 4, and 5. The appropriate as- 
sumptions and modified equations will be presented in the appropriate chapters. 



CHAPTER 3 


CIRCULAR RESTRICTED PROBLEM OF THREE BODIES 


For the planar circular restricted problem of three bodies , the 
motion of the massless particle takes place in the plane of motion of the 
two primaries and the primaries are each in circular orbits about their 
center of mass. This leads to the following reduced equations. 

3.1 Equations for the Circular Restricted Problem of Three Bodies 


The equation of motion is the same as Equation (2.1), and 


R 2/l = 6 x 


R 2/l = e y 


* 2/1 = 1 and r 2 /! ~ 0 ( 3,1) 


which implies that 


*1 * • S 2 = (1 - P): x 


(3.2) 


The out-of-plane component z = 0 , and 

a) = e , a) = 0 

z 


(3.3) 


The Jacobian function, J , [Equation (2.16)] remains the same, but the time 
rate of change of the Jacobian function, J [Equation (2.17)], is zero. 

The direction of the perturbing acceleration becomes 


n 

P 



, : yrp 

” / 3 
+ r 2 /r 2 


(3.4) 


for motion relative to primary one. For motion relative to primary two, the 
direction of the perturbing acceleration is 


17 
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n 

P 


(-e + r /r ) 

x 11 


-e + r /r 
x 11 


(3.5) 


The interval of propagation for the independent variable r_^ 

* 

(i = 1,2) is the same as Equation (2.25). Since the Jacobian function 

is a constant, the trapezoidal integration [Equation (2.26)] is not used. 

The linearized equations for 6r and 6r are the same as Equations (2.35a), 
and (2.35b) with the exception that 

- - 2 _ _ 3 _ 3 

B = e x r - r + (1 - y)r /r + yr 0 /r 0 (3.6) 

z 112 2 


3.2 Velocity Correction Direction n 

The velocity is corrected in the time-averaged direction of the per 
turbing acceleration over the propagation interval. Since the independent 
variable is not time but is the position vector magnitude, this direction is 
approximated. 

Figure 3.1 shows the variation of mean anomaly with the position 
vector magnitude for elliptic and hyperbolic conic orbits. Since the change 
in mean anomaly M is proportional to the change in time, these curves can 
be used to approximately determine the fraction of Ar corresponding to 
At/2 . 

Except near perifocus , the slopes of the curves shown in Figure 
3.1 increase with increasing r . However, this increase is less for hyper- 
bolic orbits than for elliptic orbits (see Figure 3.1). The average direc- 
tion of the perturbing acceleration is approximated by choosing the positions 
of the non-primary force center and the massless particle to the following: 
on the Earth side of the mean surface of influence, their positions at the 
correction point are used; on the Moon side, the position of the Earth 2 At/3 
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before the correction point and the position of the massless body 2Ar/3 
before the correction point are used. 
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Figure 3.1 The Variation of Mean Anomaly and Position Vector Magnitude 
3.3 Position Correction Direction n 

c 


Two ways of calculating n have been used. The first method uses 
the polar coordinates referenced to e^ which are used to describe r^ 9 

r 2 , and r> 2 (see Figure 3.2). In terms of the base vectors asso- 
ciated with these coordinate systems, the expression for n c is 

n = e cos d> ■ + e sin 6 (3. 7) 

era 



l 
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where cf> is an angle measured positive counter-clock-wise from e^ . Thi: 
method of choosing n c is called the "refined method" , because <f> can be 
varied to give accurate results [see Table 1]. The variation of <J> with 
perilune altitude is shown in Figure 3.3. 



Figure 3.3 The Variation of <j> With Perilune Altitude for the 
Circular Restricted Problem of Three Bodies 


3.4 Numerical Results for "Refined" Choice of n c 


In the non-dimensional system, the reference quantities are the 


following : 


reference mass 


reference length 


sum of the dimensional masses of the two pri 
maries 

dimensional distance between the two primaries 
[60.2684 Earth radii (er) for the circular 
restricted problem of three bodies 

o 1/2 

[(ref. length) /G(ref. mass)] 

104.21989489 hrs for the circular restricted 
problem of three bodies 


reference time 
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where G is the universal gravitational constant. The reference length was 
chosen as the average distance between the Earth and the Moon and the total 
mass of the Earth and the Moon. The masses of the Earth and Moon produced a 
mass-ratio parameter [Equation (2.2)] of y = .012150446995297 . 

The only fixed initial condition is 

r 1 = .0173014 ( = 147 n. mi altitude) (3.8) 


The initial value of the angle (see Figure 3.2) and the velocity com- 

• • 

ponents r^ and r^0^ are varied in order to attain different perilune 
altitudes. On the Earth side of the mean surface of influence 


Ar = .4977476 ( = 30 er ) 

o 

Ar- = .01659244 ( = 1.0 er) (3.9) 

r = .0173014 ( = 147 n. mi) 

o 

(r^ varies as the trajectory changes). 

For the Moon side of the mean surface of influence. 


Ar = -.01659244 ( = -1.0 er) 

o 

Ar f = -.03318488 ( = -2:0 er) 

r = .1659244 ( = 10.0 er) 

o 

r f = .00481180 ( = .29 er) 


(3.10) 


The values of , a 2 3 V 2^2 anc ^ " t ^ Tne P er il une are compared with the 

integrated results, obtained with a Fehlberg^^ Runge-Kutta RK 7(8), and 
with the patched-conic values. All runs were made on the CDC 6600 digital 
computer at the Computation Center of The University of Texas at Austin. 

Execution times for this method are 0.32 seconds per run compared 
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to 2.0 seconds for the integrated trajectory and 0.12 seconds for the patched 
conic method. 
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Table 1. Numerical Results at Perilune for "Refined” Choice of 
for the Circular Restricted Problem of Three Bodies 


A. 

73n. mi perilune altitude (<J> = 

47.164°) 




^2 

«2 (rad) 

r 2®2 

i) 

Integrated results 

.0048727 

.00476 

-2.47678 

2) 

"Jacobi" method 

.0048760 

.00903 

-2.47603 


Difference between 

.066% 

.00427 

+ .031% 


(2) and (1) 

( .67n. mi) 

( .245°) 

( +2 . 55fps ) 

3) 

Patched Conic Method .0057996 

.03138 

-2.28714 


Difference between 

19.0% 

.02662 

+7.656% 


(3) and (1) 

(192. 37n. mi) 

(1.527°) 

( +637 . 5fps ) 

B. 

501n. mi. perilune 

altitude (<j> = 

-50.658) 


1) 

Integrated 

.0069359 

.000155 

-2.12414 

2) 

"Jacobi" 

.0069360 

.000292 

-2.12400 


Difference 

.001% 

.000137 

+ .006% 



( . 0175n . mi) 

( .008°) 

(+ .459fps ) 

3) 

Patched Conic 

.00839063 

.032806 

-1.95088 


Difference 

20.97% 

.03566 

+8.156% 



(302n. mi) 

(2.042°) 

(+582fps) 

C. 

1009n. mi perilune 

altitude (<j> = 

-51.487° ) 






V* 

1) 

Integrated 

.0093806 

.0000007 

-1.872641 

2) 

"Jacobi" 

.0093802 

-.0019233 

-1.872487 


Difference 

-.005% 

-.00193 

+ .008% 



(-.0919n. mi) 

(-.111°) 

(+ . 52fps ) 

3) 

Patched Conic 

.01048787 

-.006722 

-1.764465 


Difference 

11 . 80% 

-.006730 

+5.78% 



(229.8n. mi) 

(-.386°) 

(+363. 6fps ) 


e 


n 

c 


time (hr) 

68.703 

68.641 

-.062 

69.773 

1.069 


73.182 

73.119 

-.063 

74.425 

1.244 


77.165 

77.107 

-.058 

78.701 

1.535 



-63.039) 




D. 20Q0n. mi perilune altitude ( <f> = 



T _2 

a 2 (rad) 

r 2®2 

time (h: 

1) Integrated 

.01415734 

.0000012 

-1.592244 

83.116 

2) "Jacobi" 

.01416068 

.0011052 

-1.591828 

83.007 

Difference 

.021% 

.0011039 

+ .026% 

-.108 


( .608n. mi) 

( .063°) 

(+1.40fps) 


3) Patched Conic 

.01608699 

-.000542 

-1.489750 

84.881 

Difference 

13.63% 

-.00054340 

+6.437% 

1.766 


(400. 5n. mi) 

(-.031°) 

(+344. 5fps ) 

1.766 

E. 2995n. mi perilune altitude (<\> 

= -63.395) 



1) Integrated 

.0189500 

.0000002 

-1.431530 

87.755 

2) "Jacobi 11 

.0189542 

- . 0011399 

-1.431254 

87.654 

Di fference 

.022% 

-.0011401 

+ .0193% 

-.101 


(.833n. mi) 

(-.065°) 

(+ .929fps) 


3) Patched Conic 

.0261182 

-.0294720 

-1.347823 

89.894 

Difference 

8.77% 

-.0294721 

+5.85% 

2.140 


(-1.690°) (+281 . 4fps ) 


(344.9n. mi) 
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3.5 Numerical Results for the "straight-forward" Choice of 

For the "straight-forward" choice of is 

n = an (3.11) 

c p 

where a = +1 for the Earth side of the mean surface of influence , and 

a = -1 for the Moon ^ide of the mean surface of influence. This choice 

of will be discussed in Chapter 6. The only other difference between 

the two choices of n is that on the Moon side of the surface of influ- 

c 

ence the position of primary one is At before the correction point instead 
2 

of — At and the position of the massless particle is Ar instead of 

O 

2 

— Ar . This, is due to the fact that for the Earth side of the mean surface 
of influence the direction of the perturbing acceleration [Equation 

(3.4)] is almost in the direction of the perturbing body. But for the Moon 
side of the mean surface of influence, the direction of the perturbing ac- 
celeration n^ [Equation (3.5)] leads the perturbing body by a significant 

amount. Therefore, it is more reasonable to use the position of the pri- 

2 

mary one at At instead of — At before the correction point. 

vj 

The numerical results for the straight-forward choice of n are 

c 

presented in Table 2. As shown in Table 2 the errors between the "straight- 
forward" choice of n^ and the integrated results are all less than 10 n. mi 

and 10 fps. Execution times for the choice of n are about the same as for 

c 

the "refined" choice of n^ . Therefore, it would seem logical to use the 
"straight forward" choice of n^ as it does not change with perilune altitude. 
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Table 2. Numerical 

Results? at Perilune for the 

"Straight-Forward" 

Choice oF n 


For the 

Circular Restricted Problem 

of Three Bodies 


A. 

73n. mi Perilune 

altitude 






f2 

<*2 (rad) 

r 2®2 

Time (hr) 

i) 

Integrated 

.0048727 

.004760 

-2.47678 

68.703 

2) 

" Jacobi" 

.0048615 

.004150 

-2.47902 

68.635 


Difference 

-.231% 

-.000615 

.091% 

-.068 



(-2. 33n.mi) 

(-.352°) 

( 7 . 54fps ) 


3) 

Patched Conic 

.00579909 

.031340 

-2.28722 

69.773 


Difference 

19.0% 

.026577 

- 7.65% 

1.069 



(I92.26n. mi) 

(1.537°) 

( 637 . 2fps ) 


B. 

50 In, mi perilune altitude 




1) 

Integrated 

.0069359 

.000155 

-2.12414 

73.182 

2) 

"Jacobi" 

.00693658 

.000668 

-2.12396 

73.108 


Difference 

.0092% 

.000514 

- .0086% 

-.073 



(.133n. mi) 

( .029°) 

(-.612fps) 


3) 

Patched Conic 

.0083900 

.03272 

-1.95101 

74.425 


Difference 

20.95% . 

.032566 

-8.15% 

1.244 



(301.6n. mi) 

(1.866°) 

(-581. 9fps ) 


C. 

1009n. mi perilune altitude 




1) 

Integrated 

.0093806 

.0000007 

-1. 872641 

77.165 

2) 

"Jacobi" 

.0094122 

.0009836 

-1.870040 

77.102 


Difference 

.337% 

.0009829 

- .139% 

-.0638 



(6.56n. mi) 

( .029°) 

( -8. 74fps ) 


3) 

Patched Conic 

.0111126 

.0221121 

-1.73254 

78.604 


Difference 

18.46% 

.022105 

-7.481% 

1.439 


(359. 5n. mi) (1.268°) (-470.9fps) 
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D. 2000n. mi perilune altitude 




12 

<*2 (rad) 

r 0 
2 2 

Time (hr) 

1) 

Integrated 

.0141573 

.0000001 

-1.592244 

83.116 

2) 

"Jacobi" 

.0141342 

-.0010344 

-1.592638 

83.050 


Difference 

-.163% 

-.0010349 

.025% 

-.066 



(-4.80n. mi) 

( .059°) 

( 1 . 32fps ) 


3) 

Patched Conic 

.0161907 

.0028447 

-1.487178 

84.867 


Difference 

14.4% 

.0028441 

-6.599% 

1.752 



(422. In. mi) 

( .163°) 

( -353 . 2fps ) 


E. 

2995n. mi perilune al'itude 




1) 

Integrated 

.0189500 

.0000002 

-1.431530 

87.755 

2) 

"Jacobi" 

.0189249 

-.0010122 

-1.431878 

87.771 


Difference 

-.132% 

-.0010124 

.024% 

-.043 



(-5.19n. mi) 

(-.058°) 

(1.167fps ) 


3) 

Patched Conic 

.0210521 

-.0147625 

-1.343576 

89.782 


Difference 

11.09% 

-.0147627 

-6.144% 

2.027 



(436.3n. mi) 

(-.845°) 

(-295. 6fps) 




CHAPTER 4 


ELLIPTIC RESTRICTED PROBLEM OF THREE BODIES 

For the planar elliptic restricted problem of three bodies the 
motion is still in the (e^e^) -plane. The only difference between the 
models for the circular restricted problem of three bodies and the ellip- 
tic restricted problem of three bodies is that the primaries are each in 
elliptic orbits about the CM. The choice of the eccentricity of the ellip- 
tic orbits is .0549, which is the average eccentricity of the Moon with 

[211 

respect to the Earth . The inclusion of the eccentricity of the primar- 
ies leads to the following equations. 

4.1 Equations for the Elliptic Restricted Problem of Three Bodies 

The equations of motion is the same as Equation (2.1), and 


R 2/l " R 2/l e x 


*2/1 * R 2/A + R 2/l*y <4 ' 1) 


R 2/l = 1 - e 2/l ° OS E 2/l 


R 2/l ” e 2/l Sln E 2/l /R 2/l 


where e 


2/1 


J 2/l 


= the eccentricity of primary 2 relative to primary 1 
= . 0549 

= the eccentric anomaly of primary 2 relative to primary 1, 
which is obtained from solution of Kepler T s Equation 

0 5 e 2/i : 2 * 

= angular rate of primary 2 relative to primary 1 


_ % _ / -| 2 \ 1/2 . 2 

f 2/l " 1 2/P /R 2/1 


h = - wR 2/i 5 x 


R 2 “ ^ y ^ R 2/l 6 x 


a) = cue 


a) = aie 


(4.3) 

(4.4) 

(4.5) 
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where 


-2R 


= ¥ 


2/1 


2/1 £ 
R 2/l 2/1 


= — 2oj 


2/1 

2/1 


(4.6) 


The Jacobian function J [Equation (2,16)] remains the same, but 
the time rate of change of the Jacobian function, J , [Equation (2.17)] 
becomes 


where 


1 _ 2 
J = -5 • h + y(l - y ^ 2/1 



h = r x r 


o 

r x r + oor 


2 


(4.7) 


(4.8) 


The term r(w • r) = 0 since u> is perpendicular to r for the planar 

model. The direction of the perturbing acceleration becomes 


n 

P 


- (e /R 
x 

e /R 
x 


2 

2/1 

2 

2/1 


+ r 2 /r> 2 ) 


+ V r 2 


(4.9) 


for motion relative to primary 1. For motion relative to primary 2, the di- 
rection of the perturbing acceleration is 


n 

P 


• ( - 1 ? i /r ! ) 

- ; x /R 2/1 + V r l 


(4.10) 


The interval of motion for r is the same as Equation (2.25). Since the 
Jacobian function is not constant (j f 0) the trapezoidal integration 
[Equation (2.26)] is used to predict the proper value for the Jacobian 
function at the end of the interval of motion. The linearized equations 
for 6r and 6r are. the same as Equations (2.35a) and (2.35b) except that 
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o 9 _ r> v 

B = {(w x r) - a) r + (1 - y) — + y — ) (4.11) 

r l r 2 

The velocity and position correction directions are the same as described in 
Sections (3.2) and (3.3). The variation of <j> [Equation (3.7)] with perilune 
altitude for the "refined" choice of n^ shown in Figure 4.1 



PERILUNE ALTITUDE ( M. Ml.) 

Figure 4.1 The Variation of <j> With Perilune Altitude for the 
Elliptic Restricted Problem of Three Bodies 

For the "straight-forward" method n^ was calculated as presented in Equa- 
tion (3.11). 

4.2 Numerical Results for "Refined" and "Straight-Forward" Choice of . 

The non-dimensional quantities for the elliptic restricted problem 
of three bodies are the same as those described in Section (3.4) with the 
exception that the reference length is the dimensional semi-major axis of 
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primary 2 relative to primary 1 ( = 60.2684 er). The remaining initial con 
ditions are as presented in Section (3.8). 

The values of r^ , , V 2®2 anci taTTle at P er i lune for the 

"refined 11 and "straight-forward" methods are compared with the numerically 
integrated results and with the patched conic and are presented in Table 3. 

Execution times for this problem are of the same order of magni- 
tude as for the circular restricted problem of three bodies . That is , the 
Jacobian method takes about 0.3 seconds as compared with 2.0 seconds for 
the integrated trajectory and 0.13 seconds for the patched conic method. 
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Table 3. 
Jacobi 


Jacobi - 


A. 


1) Integrated Results 

2) Jacobi - 1(* = 52.2°) 
Difference between 
(1) and (2) 

3) Jacobi - 2 
Difference between 
(a) and (3) 

4) Patched Conic 
Difference between 
(1) and (4) 


B. 


1) Integrated 

2) Jacobi - l(<j> = -49.8°) 
Difference 


3) Jacobi - 2 
Difference 


4) Patched Conic 

Difference 


Numerical Results at Perilune For 
- 1 "Refined” Choice of n 

c 

and 


2 "Straight -Forward" Choice of n^ 


59. 2n. mi Perilune 

Altitude 



^2 

a 2 (rad) 

r 2 4 5 2 

time (hr) 

.00480695 

-.15507493 

-2.463639 

69.084 

.00480630 
-.013% 
(-.134n. mi) 

-.15516424 

-.000089 

(-.005°) 

-2.463483 

.006% 

( . 525fps ) 

68.985 

-.099 

.00480564 

.028% 

( - . 275n . mi) 

-.15456276 
.00050095 
( .029°) 

-2.461228 
(- .097%) 

( -8 . 08fps ) 

69.006 

-.078 

.00527354 

9.71% 

-.15991492 

-.00483401 

-2.340174 

5.01% 

70.458 

1.374 


(96. 85n . mi) (.277°) (451.01fps) 

780. 5n. mi Perilune Altitude 


.00828199 

-.00361027 

-1.9844195 

69.951 

.00828722 

.063% 

( 1 . 0 8 5n . mi) 

-.00746075 

-.00385048 

(-.222°) 

-1.9845066 
- .004% 

(- . 293fps ) 

69.846 

-.105 

.00828288 

.011% 

( . 184n. mi) 

.00589730 
.00870000 
( . 500°) 

-1.9846983 
- .014% 

(- .937fps ) 

69.829 

-.122 

.00876841 

5.87% 

(100.96n. mi) 

-.02646631 

-.02285604 

(1.310°) 

-1.8997533 

4.27% 

( 284. 59fps ) 

71.458 

1.507 
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C. 1293.3n. mi Perilune Altitude 




V _2_ 

a 2 (rad) 

r 0 
2 2 

1) 

Integrated 

.01075232 

.06710700 

-1.8047416 

2) 

Jacobi — 1( (J> = 

-58.0°) .01075021 

.07082385 

-1.804806 


Difference 

-.020% 

.00371684 

- .004% 



(-.438n. mi) 

( .212°) 

(- . 219fps ) 

3) 

Jacobi - 2 

.01071372 

.07030673 

-1 . 8069572 


Difference 

(-.359%) 

.003199731 

. .122% 



(-8.01n. mi) 

( .183°) 

( 7.45cps) 

4) 

Patched Conic 

.01144598 

.04941106 

-1.7204624 


Difference 

6.45% 

-.01769594 

4.67% 



(I43.98n. mi) 

(-.995%) 

(283 . 29fps ) 



D. 1726. 9n. mi Perilune Altitude 


1) 

Integrated 

.01284166 

.11309066 

-1.6982265 

2) 

Jacobi — 1( <J> == 

-57.9°) .01284165' 

.11658178 

-1 . 69 83320 


Difference 

-.000% 

.00349172 

- .006% 



(-.001n. mi) 

( .200° ) 

( . 355fps ) 

3) 

Jacobi - 2 

.01283529 

.12227265 

-1.6984975 


Difference 

-.050% 

.00918228 

- .016 



(-1.32n. mi) 

( .526°) 

(-. . 911fps) 

4) 

Patched Conic 

.01371698 

.09959157 

-1.6148834 


Difference 

6.82% 

- . 01349 849 

4.91% 



(181. 68n. mi) 

(-.773°) 

(280.15fps) 



E. '3000. 3n. mi Perilune Altitude 


1) 

Integrated 

.01897652 

.20480733 

-1.5023960 

2) 

Jacobi - 1(4> = 

-55,8°) .01897642 

.20418418 

-1.5023960 


Di fference 

.016% 

-.00062315 

.020% 



(.648n. mi) 

(-.036°) 

(1. 03fps ) 

3) 

Jacobi - 2 

.01901113 

.21273488 

-1.4997524 


Difference 

.182% 

.00792754 

.176% 



(7.l9n. mi) 

( .046°) 

( 8 . 89fps ) 

4) 

Patched Conic 

.02022279 

.19094782 

-1.4176655 


Difference 

6.57% 

-.01385952 

5.64% 



(258.68n. mi) 

(-.795°) 

( 284 . 81fps ) 


time (hr) 


70.518 

70.407 

-.111 


70.406 

-.112 


72.039 

1.521 


70.978 

70.846 

-.132 


70.859 

-.119 


72.509 

1.531 


72.271 

72.135 

-.136 


72.161 

-.110 


73.945 

1.674 



34 


For the "refined” method the errors between the Jacobian method 
and the numerically integrated results are less than 1.1 n. mi and 1.1 fps. 
For the "straight-forward" method the errors at perilune are all less than 
10 n. mi and 10 fps. 



CHAPTER 5 


EPHEMERAL RESTRICTED PROBLEM OF THREE BODIES 
For the ephemeral restricted problem of three bodies the motion 
of the massless particle is not constrained to be planar and motion of the 
primaries is not two-body motion about the center of mass. The position and 

|* 0 “j 

velocity of the primaries are taken from the JPL ephemeris tape at each 
point of interest. This leads to the following equations. 

5.1 Equations for the Ephemeral Restricted Problem of Three Bodies 

The non-dimensional equation of motion for the massless particle 
is the same as Equation (2,1) 9 and 


R 


2/1 


R 2/l®x 


R 2/l ~ R 2/l 6 x + R 2/l f 2/l S y 


(5.1) 


where ^ 2/1 an< 3 ^2/1 are " ta ^ :en ^ vom ^h e ephemeris tape. The angular velo- 
city and acceleration of the primaries relative to the CM is determined by 


co “ we 


f 2/i e z 


R 2/l X R 2/l 


2/1 


(5.2) 


co = toe 


f 2/l S z 


R 2/l X R 2/l 2R 2/1 ( - 

,2 „3 V 2/1 2/1) 


(5.3) 


R' 


2/1 


R 


2/1 


Since the acceleration ^ 2/1 no ‘ t: gi ven » and is not easily calculated, it 

is assumed that at each instant in time that 

R, 


R. 


2/1 


2/1 


2/1 


(5.4) 


therefore , 
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The out-of-plane components, z and z , are not necessarily zero for the 
ephemeral restricted problem of three bodies. 

The Jacobian function J and the time rate of change of the 
Jacobian function J remain the same as Equation (2. .16) and (2.17) respec- 
tively. The direction of the perturbing acceleration for motion relative 
to primary one and primary two remains the same as Equations (2.23) and 
(2.24) respectively. 

t The interval of motion for the independent variable r is the 
same as Equation (2.25). The estimate value of the Jacobian function, J , 
at the end of an interval of motion is predicted using Equation (2.26). The 
linearized equations for 6r and 6r are the same as Equations (2.35a) and 
(2.35b), where A and B are defined in Equations (2.36) and (2.37). 

For the circular and elliptic restricted problems of three bodies, 
the results were presented for a "refined" choice of n^ , and for a "straight 
forward" choice of n . Since the results for the "straight-forward" choice 
of n^ were quite satisfactory, it was decided to only present the results 
for the "straight-forward" choice of n^ for the ephemeral restricted prob- 
lem of three bodies. 

5.2 Velocity Correction Direction n 

( : : E 

The velocity is corrected in the time-averaged direction of the 
perturbing acceleration over the propagation interval. Since the independent 
variable is not time but is the position vector magnitude, this direction is 
approximated by choosing the positions of the non-primary force center and 
the massless particle to be the following: on the Earth side of the mean 
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surface of influence, their positions at the correction are used; on the Moon 
side, the position of the Earth At before the correction point and the pos- 
ition of the massless body Ar before the correction point are used. The 

unit vector n is the direction from the massless body in the direction of 
P 

the perturbing acceleration. The reason for this choice of positions of the 
bodies for computing n^ rather than the choice used in Section (3.2) is not 
entirely obvious. For motion relative to primary one the direction of the 
perturbing acceleration is in the approximate direction of the perturbing 
body. Thus, it is logical to use the position of the massless particle and 
perturbing body at the time and position of the correction. For motion rela- 
tive to primary two the direction of the perturbing acceleration leads the 
direction of the perturbing body. This coupled with the "straight-forward" 

choice of n lead to the choice of the position of massless particle and 
c 

perturbing body described above. 

5.3 Position Correction Direction n c 

For the ephemeral restricted problem of three bodies only one choice 

of n has been used. This is the "straight-forward" choice as described 
c 

in Section (3.5). One advantage to this choice of ii c is that it is not de- 
pendent upon the perilune altitude. An explanation of the reason the "straight 

forward" choice of n works is presented in Chapter 6. 

c 

5 .4 Numerical Results for the Ephemeral Restricted Problem of Three Bodies 

For the ephemeral restricted problem of three bodies the non-dimen- 
sionalizing reference quantities are the following: 

reference mass = sum of dimensional masses of the two primaries 
reference length = semi -major axis of primary two relative to- 
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primary one at the epoch time (April 11, 1970, 


21 h 53 mi " 48 . 966 SSC ) 


reference time 


= 59.6789353907 er 

o 1/2 

= [(ref. length) /G(ref. mass)] 


= 102.694503141307 hrs 


[ 6 ] 

The mass of the Earth and Moon taken from the ephemeris tape 1 produced a 
mass ratio parameter [Equation (2.2)] of y = .01215052064981 . 

The only fixed initial condition is 

r = .01754507908 ( = 162. lOn. mi altitude) 

The remaining four initial conditions are varied to attain different peri- 
lune altitudes. The coordinate system used for the ephemeral restricted 
problem of three, bodies is a cylindrical system measured from the instantan- 
eous Earth-Moon plane and . from the Earth-Moon line. See Figure 5.1. 
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The following equations are true for the cylindrical coordinate system. 


r i = R i e Rl + z i e 5l 


(5.7) 


r l = R l e R 1 + R l <S l e 0 1 * Z l e 5l 


(5. a) 


where should not be confused with measured from the CM to primary 

one. On the Earth side of the mean surface of influence 


Ar 

.4977476 

( - 30 

er ) 

0 


Ar f = 

.01659244 

( = 1.0 

er ) 

r = 

.017545079 

( = 162 

.In. 


(r_p varies as the trajectory changes). 

For the Moon side of the mean surface of influence the values of 
Ar^ , Ar^. , and r^_ are the same as those presented in Section (3.4). 

The values of ot^ 9 ^2^2 5 an ^ I r 2 1 an d a *t perilune 

are compared with the numerically integrated results obtained with the RK7(8) 
and with the patched conic and are presented in Table 4. The values for 

5 and are not presented but compare well with the integrated results 

Execution times for this method are 0.47 seconds per run compared 
to 7.0 seconds for the integrated trajectory and 0.22 seconds for the patched 
conic method. The reason for the significant increase in execution time for 
the Jacobi method and the numerical integration is the extra time required 
to call the ephemeris tape at each step. The patched conic method only has 
to call the ephemeris tape three times per case. 



Table 4. Numerical Results at Perilune for "Straight-Forward" Choice of n c for the 

Ephemeral Restricted Problem of Three Bodies 
A. 69. 2n. mi Perilune Altitude 




a 2 (rad) 

R 0 
2 2 

l? 2 l 

CM 

I u 

Time (hr) 

1) 

Integrated 

.0025196 

-2.4383941 

.00490270 

2.4526046 

73.185 

2) 

Jacobi - 2 

.0006733 

-2.4368515 

.00490240 

2.4510451 

73.304 


Difference between 

-.0018463 

- .063% 

-.007%. 

- .064% 

1. 199 


(1) and (2) 

(-.106°) 

(5.21fps) 

(-.07n. mi) 

( -5 . 27fps ) 


3) 

Patched Conic 

.0508569 

-2.2093458 

.00614394 

2.2192867 

74.305 


Difference between 

.0483373 

-9.39% 

25.32% 

-9.51% 

+1.120 


Cl) and (2) 

(2.769°) 

( 773. 71fps) 

(255. lln. mi) 

(-788.13fps) 

+1.120 ' 

B. 

198. 5n. mi Perilune 

Altitude 





1) 

Integrated 

.03789431 

-2.3234660 

.00553223 

2.3344527 

73.312 

2) 

Jacobi - 2 

.03598942 

-2.3218661 

.00553144 

2.3328465 

73.433 


Difference 

-.00190490 

+ .437% 

-.014% 

-.068% 

+ .121 



(-.109°) 

( .291fps) 

(-.164n. mi) 

(-5.43fps) 


3) 

Patched Conic 

.08310205 

-2.1174711 

.006848689 

2.1253222 

74.435 


Difference 

.0452077 

8.87% 

23.80% 

-8.96% 

+1.123 



(2.590°) 

( 695 . 8fps ) 

(270. 57n. mi) 

(-706.4fps) 



-F 



C. 488. 8n. mi Peril une Altitude 





, 


• 




0^ (rad) 

R 2 9 2 

h 2 | 

|r 2 i 

Time (hr) 

1) 

Integrated 

.10475594 . 

-2.1270724 

. 00694462 

2.1337672 

73.584 

2) 

Jacobi - 2 

.10307331 

-2.1271008 

.00692813 

2.1338228 

73.707 


Difference 

-.0016826 

- .001% 

. -.237% 

.003% 

t .122 



(-.096°) 

( - . 096fps ) 

(-3.39n. mi) 

( .188fps) 


3) 

Patched Conic 

.13507698 

-1.9714833 

.00824681 

i. 9766647 

74.737 


Difference 

.0303210 

+7.315% 

18.75% 

-7.36% 

+1.153 



(1.737°) 

(525. 57fps) 

(267.6411. mi) 

(-530.68fps) 


D. 

975. 6n. mi Perilune 

Altitude 





i) 

Integrated 

.18567690 

-1.9091071 

.00931305 

1.9104482 

74.262 

2) 

Jacobi - 2 

.18454745 

-1.9072579 

.00930723 

1.9086020 

74.391 


Difference 

-.0011295 

+ .097% 

-.062% 

-.097% 

+ .129 . 



(.065°) 

(6.25fps) 

( -1 . 20n . mi) 

( -6 . 24fps ). 


3) 

Patched Conic 

.21112820 

-1.7836686 

.010864705 

1.7847919 

• 75.433 


Difference 

.0254513 

+6.57% 

16.66% 

-6.57% 

+1.171 



(1.458°) 

(423 . 72fps ) 

(318. 9 In. mi) 

(-424.46fps) 



-p 

ro 



E. 2379. 7n. mi Perilune Altitude 




« 2 (rad) 

• • 

r 

2 2 

4 2 ! 

4 2 l 

. T im e ( hr ) 

1) 

Integrated 

. 33873037 

-1.5927180 

.01614465 

1.4937280 

75.119 

2) 

Jacobi - 2 

.33598404 

-1.5908406 

.01610806 

1.5918543 

75.245 


Difference 

-.0027463 

+ .118% 

-.227% • 

+.118% 

+ .126 



(-.157°) 

( 6 . 34fps ) 

(-7.52n. mi) 

(+ 6 . 33fps ) 


3) 

Patched Conic 

. 34675503 

-1.5131718 

.017912635 

1.5140690 

76.346 


Difference 

.0080247 

+4.99% 

10.95% 

+5.00% 

+1.227 



( .459°) 

( 268 . 70fps ) 

( 363 . 37n . mi) 

(+269 .lOfps) 

-1.227 

F. 

4 5 5 0 . 9n . mi Peri lune 

Altitude 





i) 

Integrated 

.45144918 

-1.3871586 

.026708638 

1.3875241 

76.669 

2) 

Jacobi - 2 

.45023525 

-1.3840171 

.026658933 

1.3843824 

76.812 


Difference 

-.0012139 

+ .226% 

-.186% 

-.226% 

+ .143 



(-.069°) 

( 10 . 61fps) 

( -10 . 22n . mi) 

( -10 . 61fps ) 


3) 

Patched Conic 

.44829182 

-1.3217430 

.028746553 

1.3220842 

78.025 


Difference 

-.0031574 

+4.72% 

7.63% 

-4.72% 

+1.356 



(-.181°) 

( 220 . 97fps ) 

(418.85n. mi) 

( -221. 05fps ) 




Using the M straight -forward" choice of for the ephemeral restricted 
problem of three bodies, it is shown in Table 4 that Jacobi-2 correction 
method maintains the position and velocity errors to within 11 n. mi and 



CHAPTER .6 


DISCUSSION OF- CHOICE OF CORRECTION .DIRECTION 

The choice of the velocity correction . direction n is straight- 

P 

forward and does not need much explanation. But, the M refined ,T choice, and 
the '-straight-forward" choice of n^ should be discussed in greater detail. 

To do this a little background history * is in order. 

6.1 Background History 

When the Jacobian correction method was first conceived, it was be- 
lieved that if the velocity were corrected back to the proper value often 
enough there would be no need for a position correction. This method was 
first applied to the circular restricted problem of three bodies and, with 
only one scalar function,- it was thought that only one quantity could be cor- 
rected. The results of the Jacobian correction method evaluated at the term- 
inal point were a 10% to 20% improvement over the patched conic method. These 
results were compared at the point in the trajectory where the force center 
was changed from primary one to primary two. At this patch point the Jacobian 
method appeared to be much superior to the patched conic. However, the results 
deteriorated rapidly from the patch point in toward primary two. It was then 
realized that all the position errors that had been neglected along the tra- 
jectory transformed into a velocity error after switching to primary two as 
the force center. This position and velocity error at the mean surface of in- 
fluence then propagated into large errors in both position and velocity at the 
terminal point . 

6.2 Logical Steps That Led to the Choice of n^ 

This resulted in the approximation that 
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where a 


6r = aAt 

is defined in Equation (2.33). Then, 

6r = j a(At) 2 = t 6rAt 
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( 6 . 1 ) 


( 6 . 2 ) 


At this point, it was not clear whether n^ would be equal to n^ or not. 

By ignoring the difference between n^ and n^ then, 

6r = j 6rAt (6.3) 


which is the relationship necessary to obtain two quantities from one scalar 
function (the Jacobian function). The corrections Sr and Sr [Equations 
(2.35a) and (2.35b)] were easily determined from the linearized form of 
Equation (2.32). 

Since a relationship between n and n was not apparent , it 

c p 

was decided to vary the direction of through 360°. This led to the 

relationship that 


n = e cos <j> + e sin tp (6.4) 

c r oc 

as described in Section (3.3). From this assumption the results presented 
in Table 1 were generated. It was observed that, for motion relative to pri- 
mary one, the angle between n and n was always within ± 10° of zero 

p c 

and, for motion relative to primary two, the angle between n^ and n^ was 
very close to ± 180° . This seemed peculiar until the equations for n^ 
[Equation (2.23) and (2.24)] were analyzed. From Equation (2.23) 


- (fi 2/l /R 2/l * V r 2> 

! - , 3 - y 3, 

R 2/l /R 2/l + r 2 /r 2 


(6.5) 


interval a 


s r_ 


rr/“> 1^0 -pTOrMTl 



is ' clear that 


" r 2 

n I 'p rr ~ 

P 

where 

From Equation (2.24) 

n 

P 

In the interval as goes from .8 1.02 and r^ goes from .2 .02 

then 


R 


2/1 


= 1.0 


(6.7) 


( ~ 5 2/l /R 2/l * 

• R 2/l /R 2/l + /r l 


( 6 . 8 ) 


( 6 . 6 ) 


n 

P 



(6.9) 


but n^ points in front of primary one for all corrections. 

It then seemed obvious to correct the position vector in the direc- 
tion of n^ going from primary one to the mean surface of influence and, for 
the motion from the mean surface of influence into primary two, in the oppo- 
site direction to n 


n 

c 



on primary one side 
on primary two side 


( 6 . 10 ) 


The numerical results for the "straight -forward" choice of n 

c 

presented in Table 2 validated this choice of as being a good one. 

When the same choice of was applied to the elliptic restricted 

problem of three bodies and the ephemeral restricted problem of three bodies , 
the results were within acceptable limits [Tables 3 and 4] and provided 
further evidence in favor of the "straight-forward" choice. 



CHAPTER 7 


COMPARISONS, CONCLUSIONS, AND RECOMMENDATIONS 

7,1 Summary 

An approximate solution to the restricted problem of three bodies 
has been presented in the previous chapters. The theory of the Jacobian cor- 
rection method has been developed and the method has been applied to the 
circular, the elliptic, and the ephemeral restricted problems of three bodies 
with good results. Two methods of determining the position vector correc- 
tion direction (n ) have been presented. The numerical results for both 
choices of n^ have been presented for the circular and elliptic problems. 

Only the "straight-forward" choice of n^ was used for the ephemeral res- 
tricted problem of three bodies. The results have indicated that the "straight- 
forward" choice of n^ is more desirable, because it is independent of the 
perilune altitude. A qualitative comparison of methods and conclusions and 
recommendations are presented in this chapter. 


7.2 Comparison of Methods 

The purpose of this research has been to develop a method of gen- 
erating approximate Earth-Moon or interplanetary trajectories using knowledge 
of functions which are constant or slowly-varying for the exact motion. This 
method is compared with the qualitative characteristics of each of the exist- 
ing methods. 

[ 9 ] 

A comparison with the hybrid patched conic technique shows that 
the hybrid patched conic technique requires a patched conic trajectory as a 
reference trajectory. That is, a complete patched conic trajectory must be 
computed before the hybrid patched conic technique can be applied to the 
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trajectory. Although the Jacobian method use frequently corrected conic 
motion,, a complete patched conic trajectory does not have to be calculated 
before the Jacobian method can be applied. The hybrid patched conic tech- 
nique forms a curve fit to the perturbative accelerations , which is then 
analytically integrated to form corrections to the conic state vector. This 

technique works well for low perilune trajectories, but due to the nature of 

[4] 

the curve fit breaks down for high perilune trajectories . The Jacobian 
method gives good results for trajectories that have various .perilune alti- 
tudes (see Tables 1, 2, 3, and 4). 

r 4 l 

The ■ multi-conic technique has the capability of including Earth 
oblateness terms plus the effect of any number of gravitational bodies . A 
simplified computational algorithm for the multi-conic as applied to the 
restricted problem of three bodies : 

1. Advance conically in geocentric space from t to t + At . 

2. Calculate the average perturbing acceleration due to the 
geocentric motion of the moon. 

3 . Modify the state vectors at the end of step 1 by 

AV = aAt (7.1) 

AR = j aAt 2 (7.2) 

where a = average perturbing acceleration 

4. Convert the corrected geocentric state to a selenocentric state 
which is then projected back in time along the straight line 
defined by the velocity vector an amount At . This is the 
trajectory in a "no gravity field". 

5. The state is then advanced along a Moon-centered' conic an 
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amount At from the "no gravity" state. The final vector 
transformed to geocentric coordinates defines the new Earth- 
centered conic. 

6. The process is repeated using the new Earth-centered conic as 
the starting point. 

The Jacobian method consists of a conically advanced state that is 
frequently corrected which is less complicated than the multi-conic. How- 
ever, other perturbing bodies and Earth-oblateness terms have been incor- 
porated into the multi-conic method. Neither have been used in the Jacobian 
method. 

ran 

The pseudo-conic or the overlapped conic technique 1 considers 

only the effects of the two primaries in the equation of motion and these 

only within a certain distance of the Moon (a "pseudo-sphere of about 24 er 

is defined). Farther from the Moon only pure Earth conic motion is used. 

The overlapped conic technique actually takes less computer time than the 

patched conic method and the Jacobian technique requires about 2 to 3 times 
as much computational time as the patched conic. However, the overlapped 

conic technique corrects only approximately 80% of the errors from the 
patched conic method. 

A qualitative comparison between the Jacobian correction method 
and the matched asymptotic expansion technique^" 11 ,12 ,13 indicates that 
the matched asymptotic expansion technique is an analytical approximation to 
the restricted problem of three bodies , and the Jacobian method is a numerical 
approximation. The Jacobian method is much easier to formulate and program 
than the matched asymptotic expansion technique. 

7.3 Conclusions 

The investigation described in the previous chapters has been con- 
cerned with developing a method of numerically predicting, as a function of 
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time , the position and velocity of a space vehicle in the restricted problem 
of three bodies utilizing knowledge of the constant or slowly-varying func- 
tions of the motion. It has been desirable for the method to be computa- 
tionally fast and accurate as compared to numerical integrated trajectories 
and the patched conic trajectories. Numerical results indicate that the 
method provides a significant improvement over existing methods in the area 
of simplicity and is much more accurate than the patched conic method. 

Based on the results presented previously in this report , the fol- 
lowing conclusions can be drawn: 

1. The procedure for implementing the Jacobian correction method 
is very simple and straight-forward. 

2. The numerical results of the Jacobian correction method are a 
significant improvement over the patched conic and are accurate 
as compared to numerical integration. 

3. The numerical results of the Jacobian method do not deterio- 
rate for high perilune altitudes. 

4. The method of handling the slowly -varying Jacobian function 
for the elliptic and ephemeral restricted problems is adequate 
for all cases presented. 

5. The Jacobian correction method is sufficiently accurate to cal- 
culate trajectories with various perilune altitudes in the 
restricted problem of three bodies. 

6. The "straight-forward” choice of n^ , as presented in Chapters 
4, 5, and 6, yields sufficiently accurate results and does not 
depend upon the perilune altitude. 

7. The Jacobian correction method is computationally simple (and, 
therefore, fast) because it does not require a reference 
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trajectory, is not iterative, and needs no retracing. 

7.4 Recommendations for Future Study 

It is believed that further study in the following areas would be 

useful : 

1. . Implementation of oblateness effects when in the near vicinity 

of the primaries. It is recommended that the method developed 
by Penzo^^ be used as the method of incorporating the oblate- 
ness effects. 

2. Other perturbing bodies could be included. This would add 
terms to J , if the value of J at the correction point 
could predict accurately enough, a procedure similar to the 
one presented in this report could be used. 

3. Extending the method to free-return and interplanetary trajec- 
tories would .also be quite useful. 

4. It also would be useful to compare the method presented by 

r 17 1 . 

Nacozy with this method to see how similar t.he correction 


directions are. 
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